source("code/setup.R")
df <- read_csv("data/available_upon_request/whole_brain_data.csv")
df <- df |> mutate(dx=as.factor(dx))
roidf <- read_csv("data/available_upon_request/processed_ROI_data.csv") |> filter(dx==0)
roidfdx1 <- read_csv("data/available_upon_request/processed_ROI_data.csv") |> filter(dx==1)
ROI_Gcosinor <- read_rds("output/ROI_G-cosinor.rds")
ROI_Gcosinor_dx1 <- read_rds("output/ROI_cosinor_stats_BPD/ROI_G-cosinor_dx1.rds")
ROI_Gcosinor_dx01 <- read_rds("output/ROI_cosinor_stats_BPD/ROI_G-cosinor_dx01.rds")
ROI_AA <- read_rds("output/ROI_AA.rds")
ROI_AA_dx1 <- read_rds("output/ROI_cosinor_stats_BPD/ROI_AA_dx1.rds")
ROI_AA_dx01 <- read_rds("output/ROI_cosinor_stats_BPD/ROI_AA_dx01.rds")
# Merge all stats used throughout
roi_sumtab <- Reduce(function(x, y, ...) merge(x, y, all = TRUE, ...),
list(ROI_Gcosinor |> select(measure,roi,p_gc = pvalue),
ROI_AA |> select(measure,roi,p_AA = p),
ROI_Gcosinor_dx1 |> select(measure,roi,p_gc_bpd = pvalue),
ROI_Gcosinor_dx01 |> select(measure,roi,p_gc_com = pvalue),
ROI_AA_dx1 |> select(measure,roi,p_AA_bpd = p),
ROI_AA_dx01 |> select(measure,roi,p_AA_com = p)))
# Whole brain statistics
WB_Gcosinor_dx0 <- read_rds("output/WB_G-cosinor.rds")
WB_Gcosinor_dx1 <- read_rds("output/WB_cosinor_stats_BPD/WB_G-cosinor_dx1.rds")
WB_Gcosinor_dx01 <- read_rds("output/WB_cosinor_stats_BPD/WB_G-cosinor_dx01.rds")
WB_AA <- read_rds("output/WB_AA.rds")
WB_AA_dx1 <- read_rds("output/WB_cosinor_stats_BPD/WB_AA_dx1.RDS")
WB_AA_dx01 <- read_rds("output/WB_cosinor_stats_BPD/WB_AA_dx01.RDS")
wb_sumtab <- list(
WB_Gcosinor_dx0 |>
select(measure,es=rsqdm,es2=rsq,pvalue) |> mutate(test="Gcos",grp="CTRL"),
WB_Gcosinor_dx1 |>
select(measure,es=rsqdm,es2=rsq,pvalue) |> mutate(test="Gcos",grp="BPD"),
WB_Gcosinor_dx01 |>
select(measure,es=rsqdm,es2=rsq,pvalue) |> mutate(test="Gcos",grp="Combined"),
WB_AA |>
select(measure,es=chisq,pvalue=p) |> mutate(test="Fishers",grp="CTRL",es2=NA),
WB_AA_dx1 |>
select(measure,es=chisq,pvalue=p) |> mutate(test="Fishers",grp="BPD",es2=NA),
WB_AA_dx01 |>
select(measure,es=chisq,pvalue=p) |> mutate(test="Fishers",grp="Combined",es2=NA)
) |> rbindlist(use.names = T) |> rename(p=pvalue)
bsc <- rbind(read_rds("output/WB_S-cosinor.rds") |>
mutate(dx=factor(0,levels=c(0,1))),
read_rds("output/WB_cosinor_stats_BPD/WB_S-cosinor_dx1.rds"))
WB_acrodiffvar <- read_rds("output/WB_BPD_diffs/WB_acrodiffvar_test.RDS")
WB_Gcosinor <-
WB_Gcosinor_dx1 |>
select(amplitude,measure) |>
inner_join(WB_Gcosinor_dx1) |>
rename(acro_pcf = acrophase) |>
rename(p_pct = pvalue)
d2 <-
bsc |>
select(acro_sscf=acrophase,amp_sscf=amplitude,measure,p_ssct=pvalue,dx) |>
mutate(measure=factor(measure,levels = fig1_modord)) |>
arrange(measure)
plotdf <- WB_Gcosinor |> inner_join(d2,by="measure") |>
mutate(measure=factor(measure,levels = fig1_modord)) |>
arrange(measure)
# Acrophase consistency
dxres <-
bsc |>
group_by(measure,dx) |>
summarize(cohort_resultant=acro2complex(acrophase) |> mean() |> abs(),
cm = circamean(acrophase)) |>
mutate(dxlab=dx2lab(dx)) |>
rbind(bsc |>
group_by(measure) |>
summarize(cohort_resultant=acro2complex(acrophase) |> mean() |> abs(),
cm = circamean(acrophase),dxlab="Combined"))
diffacrovar_ms <- WB_acrodiffvar$measure[WB_acrodiffvar$p<0.05]
labmap <-
WB_Gcosinor_dx0 |> ungroup() |>
arrange(match(measure,fig1_modord)) |>
mutate(label=factor(measure,ordered=T,levels=fig1_modord,
labels = paste0(int2extlab[measure],
ifelse(measure %in% diffacrovar_ms,
"-Delta(abs(rho))-p<0.05","")
))) |>
select(measure,label) |> unique()
# (Partial) Source Data
fig4aSD <-
dxres |>
ungroup() |>
select(cohort_resultant,cm,dxlab,measure) |>
inner_join(
rbind(
WB_Gcosinor_dx0 |>
ungroup() |>
filter(pvalue<0.05) |>
select(measure,acrophase,acro_l,acro_u) |>
mutate(dxlab="CNTRL"),
WB_Gcosinor_dx01 |>
ungroup() |>
filter(pvalue<0.05) |>
select(measure,acrophase,acro_l,acro_u) |>
mutate(dxlab="Combined"),
WB_Gcosinor_dx1 |>
ungroup() |>
filter(pvalue<0.05) |>
select(measure,acrophase,acro_l,acro_u) |>
mutate(dxlab="BPD"))
) |>
mutate(measure=int2extlab[measure])
write_csv(fig4aSD,"output/source_data/Figure_4a_partial.csv")
Gcos_3grp_acrotick_p <-
plotdf |>
inner_join(labmap) |>
ggplot(aes(acro_sscf,1-p_ssct))+
ylab(expression(abs(rho)))+
night_rect_early()+
geom_rect(xmin=0,xmax=24,ymin=0,ymax=1,alpha=0.1,color=NA,fill="whitesmoke")+
geom_rect(xmin=0,xmax=24,ymin=-Inf,ymax=0,fill="white")+
geom_hline(yintercept = 1,color="black",size=0.2)+
geom_hline(yintercept = 0.5,color="black",size=0.2,linetype=3)+
annotate(geom="segment",x=seq(0,24,3),y=2.2,yend=0,xend=seq(0,24,3),
color="grey90",size=0.3)+
geom_text(aes(color=dx2lab(dx),y=1,
angle=-acro_sscf/24*360-90),
size=5,label="-")+
# G-cosinor dx1
geom_point(data = plotdf |> filter(p_pct<0.05) |> inner_join(labmap),
aes(acro_pcf,y=1.6,color=ifelse(p_pct<0.05,"BPD",NA)),size=1.2)+
geom_errorbar(data = plotdf |> filter(p_pct<0.05) |> inner_join(labmap),
aes(xmin=acro_u, xmax=acro_l,y=1.6,
color=ifelse(p_pct<0.05,"BPD",NA)),width=0.25)+
# G-cosinor dx0
geom_point(data = WB_Gcosinor_dx0 |> filter(pvalue<0.05) |> inner_join(labmap),
aes(acrophase,y=1.3,color=ifelse(pvalue<0.05,"CNTRL",NA)),size=1.2)+
geom_errorbar(data = WB_Gcosinor_dx0 |> filter(pvalue<0.05) |> inner_join(labmap),
aes(xmin=acro_u,x=acrophase,
xmax=acro_l,y=1.3, color=ifelse(pvalue<0.05,"CNTRL",NA)),width=0.25)+
# G-cosinor dx0dx1
geom_point(data = WB_Gcosinor_dx01 |> filter(pvalue<0.05) |> inner_join(labmap),
aes(acrophase,y=1.9,color=ifelse(pvalue<0.05,"Combined",NA)),size=1.2)+
geom_errorbar(data = WB_Gcosinor_dx01 |> filter(pvalue<0.05) |> inner_join(labmap),
aes(xmin=acro_u,x=acrophase,
xmax=acro_l,y=1.9, color=ifelse(pvalue<0.05,"Combined",NA)),width=0.25)+
geom_segment(data=dxres |> inner_join(labmap),
aes(xend=cm,x=cm,y=0,yend=cohort_resultant,color=dxlab),
size=0.5,alpha=0.8,
arrow = arrow(type = "closed",length = unit(0.8,"mm")))+
scale_color_manual(values=c("BPD"=colors$dx1,
"CNTRL"=colors$dx0,
"Combined"=colors$dx01),guide='none')+
facet_wrap(label~.,nrow=3,labeller = label_parsed)+
theme(panel.border = element_blank(),
panel.background = element_blank(),
panel.grid = element_blank(),
panel.spacing.x = unit(0,"line")) +
coord_polar()+
theme(panel.grid = element_blank())+
theme(panel.spacing = unit(0, "lines"))+
scale_y_continuous(breaks = c(0,1),limits = c(-0,2.2))+
scale_x_continuous(limits = c(0,24), breaks=seq(0,24,6))+
xlab("Acrophase (hr)")+
theme_condense()+ theme(panel.border = element_blank(),
panel.background = element_blank(),
#Legend
legend.title.align = 0,
legend.position = "right",
legend.title = element_text(size=8, vjust = 0, hjust = -3),
panel.spacing.x = unit(0.5,"line")
)
Gcos_3grp_acrotick_p
Extract stats referenced in text for “new” ROIs identified by combined (“FALSE TRUE”).
G-cosinor
roi_sumtab |>
group_by(measure) |>
mutate(class=paste(p.adjust(p_gc,"fdr")<0.05,p.adjust(p_gc_com,"fdr")<0.05)) |>
group_by(measure,class) |> count() |>
filter(class=="FALSE TRUE")
# A tibble: 3 × 3
# Groups: measure, class [3]
measure class n
<chr> <chr> <int>
1 cbf FALSE TRUE 17
2 fa FALSE TRUE 8
3 gm_md FALSE TRUE 117
AA
roi_sumtab |>
group_by(measure) |>
mutate(class=paste(p.adjust(p_AA,"fdr")<0.05,
p.adjust(p_AA_com,"fdr")<0.05)) |>
group_by(measure,class) |> count() |>
filter(class=="FALSE TRUE")
# A tibble: 4 × 3
# Groups: measure, class [4]
measure class n
<chr> <chr> <int>
1 cbf FALSE TRUE 119
2 gm_md FALSE TRUE 101
3 md FALSE TRUE 2
4 qt1 FALSE TRUE 93
Or the simple difference in #ROIs
roi_sumtab |>
group_by(measure) |>
summarize(n_gc=sum(p.adjust(p_gc,"fdr")<0.05),
n_gc_com=sum(p.adjust(p_gc_com,"fdr")<0.05),
n_AA=sum(p.adjust(p_AA,"fdr")<0.05),
n_AA_com=sum(p.adjust(p_AA_com,"fdr")<0.05)) |>
mutate(diff_gc=n_gc_com-n_gc,
diff_AA=n_AA_com-n_AA)
# A tibble: 9 × 7
measure n_gc n_gc_com n_AA n_AA_com diff_gc diff_AA
<chr> <int> <int> <int> <int> <int> <int>
1 cbf 289 265 128 246 -24 118
2 ct 0 0 0 0 0 0
3 fa 0 8 0 0 8 0
4 gm_md 53 168 36 136 115 100
5 md 3 0 1 3 -3 2
6 qt1 0 0 4 97 0 93
7 sa 0 0 0 0 0 0
8 vol 0 0 0 0 0 0
9 wm_qt1 0 0 0 0 0 0
library(ggseg)
library(ggsegGlasser)
meta <- get_glasser_spatial_meta()
predat <- roi_sumtab |>
pivot_longer(c(p_gc,p_AA,p_AA_bpd,
p_AA_com,p_gc_bpd,p_gc_com)) |>
group_by(name,measure) |>
mutate(q=p.adjust(value,"fdr")) |>
mutate(name=factor(c("p_gc"="CNTRL Gcos",
"p_AA"="CNTRL Fishers",
"p_gc_bpd"="BPD Gcos",
"p_AA_bpd"="BPD Fishers",
"p_AA_com"="Combined Fishers",
"p_gc_com"="Combined Gcos")[name],
levels=c("CNTRL Gcos","CNTRL Fishers",
"BPD Fishers","BPD Gcos",
"Combined Fishers","Combined Gcos"),
labels = c("CNTRL Gcos (n=16)","CNTRL Fishers (n=16)",
"BPD Fishers (n=8)","BPD Gcos (n=8)",
"Combined Fishers (n=24)","Combined Gcos (n=24)"))) |>
group_by(measure,name) |>
mutate(n=sum(q<0.05),ntot=n()) |>
filter(n>0) |>
mutate(value=ifelse(q<0.05,value,NA))
# Don't forget about these
predat |>
select(n,ntot,name,measure) |>
unique() |>
filter(measure %in% c("wm_qt1","fa","md"))
# A tibble: 5 × 4
# Groups: measure, name [5]
n ntot name measure
<int> <int> <fct> <chr>
1 5 46 BPD Gcos (n=8) fa
2 8 46 Combined Gcos (n=24) fa
3 3 46 CNTRL Gcos (n=16) md
4 1 46 CNTRL Fishers (n=16) md
5 3 46 Combined Fishers (n=24) md
library(ggsegGlasser)
library(ggseg)
meta <- get_glasser_spatial_meta()
meta_glass <- get_glasser_spatial_meta()
fig4bSD <-
ROI_Gcosinor_dx01 |>
filter(measure =="gm_md") |>
mutate(acrophase=ifelse(qvalue<0.05,acrophase,NA)) |>
select(measure,roi,qvalue,acrophase)
write_csv(fig4bSD |> mutate(measure=int2extlab[measure]),
"output/source_data/Figure_4b_part1.csv")
gmmdcombgc <-
inner_join(fig4bSD,meta) |>
ggplot(aes(fill=acrophase))+
facet_wrap(~int2extlab[measure],nrow=5)+
theme_condense()+
scale_fill_clock(rot=15,na.value="grey95")+
geom_brain(atlas = glasser,
position=position_brain(side ~ hemi),
color="grey",size=0.1)+
theme_void()
gmmdcombgc
# WM-FA Gcos comb
library(ggsegICBM)
library(plotly)
library(ggseg3d)
meta <- get_icbm_spatial_meta()
# Source Data
fig4bSD2 <-
ROI_Gcosinor_dx01 |>
filter(measure=="fa") |>
mutate(acrophase=ifelse(qvalue<0.05,acrophase,NA)) |>
select(measure,roi,qvalue,acrophase)
write_csv(fig4bSD2 |> mutate(measure=int2extlab[measure]),
"output/source_data/Figure_4b_part2.csv")
pdf2 <- inner_join(fig4bSD2,meta) |> icbm_preprocess()
pdf2$region %in% icbm_3d$ggseg_3d[[1]]$region
[1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[21] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[41] TRUE TRUE TRUE TRUE TRUE
rot=15
pal <- 1:24
names(pal) <- rainbow(24)[((1:24)-1+rot)%%24+1]
p3d <- pdf2 |>
select(region,acrophase,roi) |>
ggseg3d(atlas = icbm_3d,colour="acrophase",
palette = pal,na.alpha = 0.01) |>
add_glassbrain(opacity=0.05) |>
pan_camera("right lateral") |>
remove_axes()
p3d